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ABSTRACT 

Current research at Rensselaer, is generating 
fundamental engineering design techniques and concepts 
for the optimization of a gas chromatograph-mass spectro- 
meter chemical analysis system suitable for use on an 
unmanned, martian roving vehicle. Previously developed 
mathematical models of the gas chromatograph are inade- 
quate for predicting peak heights and spreading for some 
experimental conditions and chemical systems.. Consequently, 
a modification to the existing equilibrium adsorption 
model is required. The Langmuir isotherm replaces the 
linear isotherm. A olosed-form. analytical solution to 
the model is not available so the numerical technique of 
Crank-Nicolson is studied. Initial work considers the 
linear isotherm to determine the utility of the method. 
Modifications are made to the method to eliminate unnec- 
essary calculations. These modifications result in an 
overall reduction of the computation time of about 42%. 

After successful utilization of the Crank-Nicol- 
son method to the linear isotherm, interest is shifted to • 
the Langmuir isotherm which takes into account the com- 
position-dependent effects on the thermodynamic parameter, 
mRo. 

This model shows the sharp rise of the peak and 
the spreading of the tail generally observed in experimental 



data. A secant method is incorporated into the model to 
determine the Langmuir constant by matching the simulated 
peak time with the experimental output peak time. This 
secant scheme coupled with the use of an mRo corresponding 
to a dilute sample, simulates the actual data quite well. 

A problem is encountered with oscillations of some of 
the simulations. This is due to the difficulties in using 
a cubic spline interpolation scheme capable for handling 
a limited number of input data points. These oscillations 
also produce problems when the secant method is used in 
determining the Langmuir constant. Convergence, if it occurs, 
is often slow. In fact, for a system which has a large 
number of input data points, no convergence was achieved. 

An area of future research would be to find a 
better interpolation scheme to handle a large set of input 
data points. Further research could be incorporating a 
numerical scheme to vary mRo which in turn determines the 
Langmuir constant. In this ’•'ay the simulation curve is 
fit to the actual data by varying two constants, mRo and K. 
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PART 1 


INTRODUCTION AND SUMMARY 

One important phase of the unmanned martian 
roving vehicle is the search for and analysis of organic 
matter and living organisms on the Martian surface. The 
overall task is to provide design criteria and engineering 
techniques for use in optimizimg the design of a system. 

The plan is to use a gas chromatograph-mass spectrometer 
system, the gas chromatograph being used to separate the 
mixtures and the mass spectrometer to analyze the chemical 
compounds. There must be some flexibility in the system 
in order to perform a diversity of experiments, including 
analysis of the atmosphere and samples from incubation 
tests of soil and atmosphere. However, the most stringent 
design requirement for the system is that it be small 
and light enough to fit in the payload of the vehicle. 

This report describes the recent work done on developing 
a nonlinear simulation model to better represent the gas 
chromatograph. 

Gas chromatography is a way of separating a 
mixture of different chemical species by utilizing the 
mechanism of adsorption-desorption. Owing to the different 
characteristics of various chemicals, each species will 
adsorb and desorb at different rates when exposed to a 
packed bed of granular particles with or without a liquid 
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substrata. The more strongly a species is adsorbed, the 
longer its elution time will be. because of the unique 
behavior of each chemical, a multicomponent sample may 
be injected into a chromatograph and elute as a series 
of pulses each containing one of the components. 

Since gas chromatography is a complex process 
with many dynamic mechanisms continually taking place, 
mathematical models are being developed to predict system 
behavior for varying degrees of complexity. These 
predictions are then compared to actual data and if 
deviations from the experimental results occur then there 
is a flaw in the model. For an exact model which incor- 
porates all observed phenomena, there would' be no deviations; 
however, the complexity of such a model would make its 
numerical solution impractical even with present day com- 
puters. 

One must be satisfied with models describing 
only the significant phenomena affecting the behavior of 
the column. If consistent deviations occur from the actual 
data then an important mechanism has been neglected. 

The model then must be amended. 

Prior to this investigation, chromatographic 
models have been developed which have simulated success- 
fully many systems well. However, in a few instances, 
model predictions have deviated consistently from the 
actual data. For dilute samples the model simulates the 
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system well but for some chemical species deviations 
occurred as the sample was increased. It is the objective 
of this report to develop and evaluate a practical nonlinear 
adsorption model to account for composition-dependent 
effects observed in previous research. 

Two versions of the equilibrium adsorption model 
have been investigated. The first version, Equation Ilia 
of Figure 1, uses a linear isotherm. With the linear isotherm 
analytical techniques can be used to find the solution 
of the resulting second-order equation. The second ver- 
sion, Equation Illb of Figure 1, incorporates the Langmuir 
isotherm into the model as a means of accounting for the 
observed composition-dependent effects. Analytical tech- 
niques are not . available for solving the resulting non- 
linear equation so numerical methods are being used. 

Because of instabilities with an orthogonal collocat ■ . 

procedure used previously (Lavoie, 1974), the simple, 
finite difference scheme of Crank-Nicolson is being used. 

This method was first applied to the linear version and 
the results were compared with the available analytical 
results. Considerable effort was made in coding to elim- 
inate calculations involving multiplication by zero in 
order to reduce computation time. This additional coding 
reduced the computing time by about 35%. Good comparison 
betv/een the numerical , and analytical results gave confi- 
dence in the method. 



In applying the method to the nonlinear ver- 
sion, a modification was required because of the nonlinear 
term. A predictor-corrector technique was used. Initial 
work with this model used the value of the thermodynamic 
parameter, mRo, obtained by a curve fitting technique 
between the simulated output and the actual output. The 
peak times of the simulation output was compared with 
the actual peak time and there was fairly close agreement. 
However, when the Langmuir constant, K, was increased 
the peak times occurred earlier. The peaks for the simu- 
lated output were much larger than the actual peaks (about 
90% greater) . Qualitatively, as the Langmuir constant, 

K, was increased the peaks became smaller approaching 
the actual value. More importantly, the shape of the sim- 
ulated output appeared more like the actual shape for 
higher values of K. The curves rose more iharply and exhib- 
ited more tailing as K was increased. 

Something had to be done to alleviate the prob- 
lem of the early occurrence of the peaks. One possibility 
substantiated by the fact that mRo was composition-depen- 
dent was to correct tor the fact that the sample is not 
dilute. Extrapolate the value of mRo to a very dilute 
sample and use this for the larger sample. In this way 
the Langmuir constant can be varied until the peak times 
are closely matched. After these peak times are matched, 
a comparison of their peak heights can be made. A sim- 



ulation peak height greater than the actual peak height 
suggests that a smaller mRo should be attempted next. 

A simulation peak height smaller than the actual peak 
height would indicate a larger value of mRo be used. This 
procedure is in fact a two constant curve fitting tech- 
nique to obtain a simulation curve similar to the exper- 
imental data. 

When this procedure was incorporated into the 
program close agreement resulted between the simulated 

v 

chromatograph and the actual chromatograph. Especially 
for a nondilute sample, the Langmuir isotherm with this 
two constant curve fitting simulates the actual data 
quite well. 

This procedure of matching the peak times of 
the simulation with the actual data was then applied 
to the system of n-heptane on the Chromo sorb-10 2 at 175°C. 
This system's input had a large set of data points so 
that oscillations in the peak occurred. Due to these 
oscillations the program had" difficulties converging. 


PART 2 


BACKGROUND 

One area of the overall gas chromatograph systems 
investigation has been the mathematical modeling of the 
chromatograph system. Many previous investigators have 
contributed to the development of mathematical models 
(Keba and Woodrow, 1972; Meisch,1973; Sliva,1968; Taylor, 

1970; Voytus, 1969). A course has been pursued wherein 
successively more complex models have been considered. 

These models have all yielded analytical expressions from 
which a simulation chromatograph could be computed directly. 

Ccwapari-so-n. ©-£ predicted system behavies with, a ctuaX .. sy stem 

data has directed modeling efforts to consider more ade- 
quate and hence more complicated models. 

A system of three coupled partial differential 
equations was derived earlier to model the gas chromatograph 
(Woodrow, 1974) . Excessive computation time due to the 
complexity of the system of equations indicated a simpli- 
fication of the model was required. Voytus (1969) derived 
a simplification of the model which considered axial dif- 
fusion, convection, an <j equilibrium adsorption/desorption 
using the linear isotherm. Assumptions for the model, 
designated the equilibrium adsorption model, included 
negligible mass transport effects between the carrier gas 
and the adsorbent phase and no intra-particle diffusion 
effects. In its mathematical fo^m, the model is represented 
by the following second-order, linear partial differential equation 

~i \ 
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which was derived from the mass balance for the injected 
species in the fluid and solid phases of the column: 


(l+l/mRo)^y/ie = (1/Pe) d 2 y/dz 2 - dy/^z 


} D 

T s 
P . 

>• 0 


Boundary conditions: 
y(O,0) = input pulse 
lim y(z,9) = finite 

2 -* 09 

The final solution of this equation is the following 


convolution integral: 


f 9 

:/©) = J y (0,t-r)- 


y f (x)- dt 


where y(0,e) = the input pulse 

frl = the unit impulse response at the end 

of the column 

= 1 / gPe expj-Pe (?-/*)> 

2\/Trr 1 4 JF J 

and f) - (1+1/nRo) 

The Peclet number is a measure of the dispersion of the 
sample component due to the axial diffusion in the column. 
Smaller values of Pe indicate greater dispersion. The 
thermodynamic parameter, mRo, determines the elution 


time of the injected component. Larg|r values of mRo 
indicate that the species are not strongly adsorbed and 
hence they will elute at an earlier time. A study of the 
nonequilibrium adsorption model has shown that mass trans- 
port between the carrier gas and the adsorbent is too 
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fast to have any limiting effects on systems studied / 
thus far (Keba and Woodrow, 1972) . 

For many systems the equilibrium adsorption 
model with the linear isotherm predicted quite well the 
results obtained experimentally. In general, these cases 
were dilute sample sizes end particular chemical species. 
For other chemical species, the model predictions deviated 
consistently from actual data as the sample size was in- 
creased. Figure 2 compares the actual output data to the 
model simulation for n-heptane on Chromosorb-102 at 175°C 
and there is little resemblance between the two. The 
actual data is very spread with much tailing, while the 
modeled response is a relatively sharp peak with little 
tailing. 

Meisch (1973) found that the parameter, mRo, ^ 
was a significant function of composition. This was 
further substantiated by Lavoie (1974) with experiments 
with n-heptane and a variety of sample sizes. Figure 3 
shows results of those experiments. The parameter, mRo, 
decreased linearly with a decrease in sample size. 

The equilibrium adsorption model with the linear 
isotherm obviously does not represent all of the physical 
behavior that is occurring in the column. There are many 
possibilities to consider. The variation of mRo with 
composition supports the idea that a nonlinear isotherm 







Figure 2 Comparison of Actual Chromatograph and Equilibrium Adsorption 
Model {Linear Isotherm) for n-Heptane on C-102 at 175°C 
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should be used in the model. Preliminary work (Lavoie, 
1974) has already begun using orthogonal collocation to 
solve the resulting nonlinear equation but due to the 
instabilities in the results, a new different numerical 
scheme .is being proposed, Crank-Nicolson. This idea is 
investigated in this paper. 
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PART 3 


RESULTS AND DISCUSSION 

A . Formulation of the Crank -Nicolson Algorithm to tho 

Equilibrium Adsorption Model Using the Linear Isotherm 

One mathematical model that has already been 
developed is the equilibrium adsorption model* using the 
linear isotherm to describe the adsorption kinetics. This 
model considers interparticle axial diffusion, convection 
and equilibrium adsorption/dcsorption. Diffusion of the 
chemicals in the direction of the carrier gas flow in the 
interparticle region is represented by the- dimensionless 
parameter, Pe, which is determined by the system fluid- 
mechanics. Adsorption-desorption on the particle is rep- 
resented by mRo, a thermodynamic parameter peculiar to 
each species. This parameter contains an equilibrium 
constant, m, and the quantity, Ro. Ro is the ratio of 
the moles of fluid within the total bed to the moles of 
adsorptive sites within the bed. 

With the above concepts in mind, the following 
dimensionless equation results; 



With boundary conditions; 

y(O,0) = input pulse (2) 


1 For ”tHe details of the derivation see Appendix A 


m~.m * » , 
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lim y(z,9) = finite (3) 

z -* » 

and initial condition: 

y(z,0) = 0 This initial condition states 

that initially the column has 
only inert carrier gas flowing 

through it. (4) 

An alternate terminal boundary condition that 
could be used instead of the finite composition at an 
infinite distance down the column (equation (3)] is the 
following proposed by Woodrow (1974) : 

$y(l,e) = 0 for 9>0 (3a) 

3 ^ 

Use of the infinite boundary condition (equation 
(3)3 in analytical work yields a great deal of mathematical 
simplification and is consistent with the theory which has 
been developed. However, when numerical techniques must 
be applied to solve equation (1) , the terminal boundary 
condition given by equation (3) must be replaced by a 
terminal boundary condition which is both computationally 
expedient and physically meaningful. The boundary con- 
dition proposed is given by equation (3a). 

To obtain numerical solutions to partial dif- 
ferential equations, continuous variables are replaced by 
discrete variables. The relations between these discrete 
variables in the method of finite differences are finite 
difference equations which are based on Taylor series 





representations of the dependent variable. The domain of 
the independent variable that is discretized form a system 
of grid points. Figure 4 shows a grid representation for 
the transient analysis of a system. The spatial dimension, 
z, is shown as being bounded and normalized and the time 
domain, 0, is shown with no particular bound. The grid 
is fixed; i.e., spatial discretizations and time discre- 
tizations are uniform for each domain. The value of z, 
the continuous space dimension is given by: 

z — i * (Lz) 

where i refers to a particular spatial grid point and 
4z is the spacing between grid points. Similarly, the 
value of 0, the continuous time dimension is given by: 

0 = n •’(A0) 

where n refers to a particular time grid point andA© 

is the interval between grid points. 

For the parabolic system of the second order 

chromatograph system model, equation (1), the two-level 

implicit method known as the Crank -Nicolson method of 
I 

solution is probably most popular and is well documented 

i 

(Lapidus, 1962; von Rosenberg, 1969). The Crank-Nicolson 
method replaces the partial derivatives which are eval- 
uated at the position (space i, time n+M with their 
finite difference counterparts. These approximations are 



a 
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an arithmetic average of the finite difference analogs at 
the points z^, 0 n and z^, © n+ i. The approximations at 
each of these points can be viewed also as as arithmetic 
average so that the resulting analog takes on the form 
of a double average of the compositions surrounding the 
point, y i>n+/t . In this method, the following approximations 
are made for the first and second spatial derivatives 
and the first time derivative: 



These finite differences are substituted into 


the partial differential equation, giving a system of 

* 

algebraic equations having a tridiagonal system matrix. 

This system of equations is solved using the Thomas al- 
gorithm (von Rosenberg, 1969) . 

Preliminary studies were done with this procedure 


0 

II 

II 

n 


to determine the appropriate choice of the incremental 
variables, Ae and Az. Von Rosenberg (1969) has determined 
that the Crank-Nicolson method becomes unstable (oscillatory 

1 FoiT "detaiTs~oT ' de r I va t i on see Appendix B 
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about the true curve) if the coefficient of the second 
space derivative is small compared to the coefficient of 
the first space derivative. There will be no oscillations 
when 

coeff. of first space dor. x 4 1 
coeff. of* second space der7 2 

For the equilibrium adsorption model, equation (1) reduces 
to 

Pe -Az £ 1 
2 

For practical chromatography, the Peclet numbers 
are in- th-e- order of 10 , 000 . computer simulations 

for Az less than the values specified by von Rosenberg 
exhibited oscillations. For a z « 0.001 no oscillations 
were..found. All of the remaining simulations were obtained 
with this value of the incremental space dimension. 

For the Crank-Nicolson algorithm there is no 
stability criterion for the ratio of the incremental vari- 
ables. Having determined the value of Az needed, theoretically 
any time increment could be used. The effect of time 
step size was investigated. Table 1 shows the computation 
time needed for each value of Ae. For small values of 
AQ the computation time is excessive. (Note the computa- 
tion times in the table are for equal partial simulation 
runs.) For A9 smaller than 0.04 the decrease in computation 
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TABLE 1 


COMPUTATION TIME FOR VARIOUS 
DIMENSIONLESS TIME INCRExMENTS 


0 

n 


Dimensionless Time 
Increment, Q 


Computation Time for 
Equivalent Partial Run (sec.) 


U 


0.0004 

1840 

0.004 

334 

0.008 

204 

0.01 

171 

0.04 

87 

0.05 

87 
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time was almost negligible while the magnitude of the 
' differences from the values calculated for smaller d© 
increased. As a result of this study the value of 0.04 
was used for the incremental time dimension. 

B . Modification of the Crank-Nicolson Method 

When the decision was made on the incremental 
variables, modifications were necessary to further reduce 
the computation time. The modifications were made to 
the Thomas algorithm (von Rosenberg, 1969) . This algorithm 
is used to solve a tridiagona t. system of algebraic 
equations. of the form: 

a^i-l + b i^i + c iyi+l = djL (8) 

for 1 -i -R 
v/ith a^ = c R = 0 

The algorithm is as follows: 

First, compute 

h = b i " a i c i-l with " b l O) 

Ji-1 . 

and 

tfi = di - with tfj. = d^ (10) 

h b l 

The values of the compositions (y) are then computed from 
y R = tf R and y L = - c iYi+1 

h 


( 11 ) 


For early times in the simulation, the pulse 
has not had an opportunity to spread out throughout the 
column. Since the 'd' vector [equation (8)] is a linear 
combination of the compositions at the spatial points 
(i-1, i, and i+1) for the previous time step (n) , many 
of its elements are zero. The trend of thought is why 
do all the multiplications by zero when this will just 
result in zero compositions for the next time step (n+1) . 
Because values of ’d' are either very small or zero, they 
may produce very small values for an intermediate calcu- 
lation vector, of equation (10). When the ith element 
of the 'd' vector coupled with the ith element of the ' 
vector are both small, then the corr e-spondi rrg ith element 
of the composition vector is set equal to zero. 

A numerical scheme was developed to determine 
when to eliminate unnecessary calculations. Between the 
indices, II and 12, are the positions at any one time 
step in the column where the pulse is located. The com- 
positions are set equal to zero to the left of II and to 
the right of 12. This feature also eliminates many of 
the underflows previously encountered. 

With this modification the computation time 
was reduced by about 35% for an entire simulation. The 
time is saved during the early part of the simulation 
before the pulse has propagated faj down the column where 
there is no fluid in the late sections of the column. 
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Time is also conserved at the end of the simulation when 
.there is no fluid in the beginning parts of the column. 

A program has been written to solve the 
equilibrium adsorption model with the linear isotherm 
using the Crank-Nicolson finite difference scheme. Figure 
5 shows the organization of the program. After the coding 
of the equilibrium adsorption model by the Crank-Nicolson 
method was completed, results obtained from the method 
were compared with the analytical solution. Figures 6, 

7, and 8 show the actual data, the analytical solution, 
and the simulation by the Crank-Nicolson method. There 
is very little difference between the Crank-Nicolson 
s im u lat ion, and the analytical solution and these two 
methods are compared for various points in Table 2. Even- 
though these two curves are in close agreement with them- 
selves they do not model actual experimental data so the 
incorporation of the Langmuir isotherm into the equilibrium 
adsorption model was investigated. 


C. Incorporation of the Langmuir Isotherm into the 
Equilibrium Adsorption Model 


, A nonlinear isotherm of the Langmuir type may 

possibly explain the deviations mentioned previously. 

The Langmuir isotherm gives the following relation between 

\ 

the adsorbed phase and mobile gas phase mole fractions 
(x and y respectively) : 



Figure 5 Crank-Nicolson Method Program for the 

Equilibrium Adsorption Model (Linear Isotherm) 
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Figure 6 The Actual Chromatograph Data for n-IIeptane on 
Chromosorb-102 at 175 # C 
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TABLE 2 


COMPARISON OF 

THE CRANK-NICOLSON METHOD (LINEAR ISOTHERM) 

WITH 

THE CONVOLUTION (ANALYTICAL) SOLUTION 


Dimensionless Crank- 


Time 

' Nicolson 

28.3058 

0.018971 

29.0047 

0.020664 

30.1696 

0.01646 

33.8971 

0.0042437 

35.5555 

0.0010353 


Convolved 

(Analytic) 

Relative 

Error 

0.0180004 

5.38% 

0.020698 

0.16% 

0.016995 

3.15% 

0.0042442 

0.012% 

0. 001(1643 

2.73% 



II 

II 

0 

0 
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x 


, *y . 

1 + Ky 


This equation is compared to the linear isotherm in Figure 
9. For small concentrations the Langmuir isotherm has 
nearly linear behavior. For large concentrations, the 
amount that can be adsorbed becomes proportionately less 
and finally approaches a maximum value. At high concen- 
trations the adsorbent sites become saturated so that an 
increase in the gas phase concentration has little effect 
on the adsorbent phase concentration. Larger values of 
K make the isotherm nonlinear at lower mobile gas phase 
concentrations. 

The nonlinearity would have the following effect 
on the behavior of the column. Large samples entering 
the column would be adsorbed proportionately less in zones 
of high concentration than in low concentration. As a 
result large samples would begin to elute sooner than 
small samples and would be more spread. The quantitative 
effects on elution time, which is related to mRo and 
spreading can be found only by solving the model. 

This model takes on the form* 

ft + _ fi ~U 2 y - h. < 12 > 

l mRo(l+Ky) 2 J 50 “ (Pej^2 6 z 


* See Appendix A for derivation 


Gas Phase Concentration 


Figure 9 


Comparison of the Linear and Langmuir 
Isotherms 
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Comparing this equation to the equilibrium adsorption 
model [equation (1)], the only difference is the term, 
(l+Ky)^ multiplying mRo. In effect, the parameter, mRo, 
becomes a variable dependent on composition, y. This 
dependence becomes stronger for increasing K. The aver- 
age value of mRo(l+Ky) 2 increases with sample concentration. 
This is the observed behavior in experimental data. 

The equation must be solved numerically because 
of the nonlinear coefficient. Based on prior research 
(Lavoie, 1974), the Crank-Nicolson method will be used 
instead of orthogonal collocation since orthogonal collo- 
cation produced instabilities. 

In applying the method to the nonlinear equation, 
a modification is required because of the nonlinear term. 

A predictor-corrector technique is used. Initially the 
partial derivatives are replaced by the finite differences 
as in the linear version. The composition-dependent 
coefficient of the time derivative is treated by a method 
developed by Douglas (1958) . The composition in the non- 
linear term has to be evaluated at the position (space i, 
time n+*/J . This value, yj.,n+v * s approximated by its 
value at the previous time step, y^ n , and a correction 
factor, (A e/2) (<*y/*e) 4 n where © is time. This derivative 
is approximated by its value from the differential equation 
in which the derivatives and compositions are evaluated 
at the previous time step. The mathematical equations 
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for this approximation of y. are: 

1 , n+ 'it 


*!,»+** y 

where 
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i,n 



1 + mRoTI+Ky i#n ) 2 



Pe 
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(13) 

(14) 


This approximation for the nonlinear term is substituted 
into the differential equation and the resulting system 
of algebraic equations is solved by the Thomas algorithm. 
The modification incorporated into the Thomas algorithm 
for the linear isotherm is again used shorter computation 
times resulting. After the first determination of the 
compositions in the column, another iteration has to be 
done. During these subsequent iterations a new estimate 
is made for the nonlinear term. This is done by averaging 
the compositions. 


y i # n+Yt _ y i,n + y i,n+l 

m* , - u ' ‘ "■ 111 J" -r *ti i_ i -|l u _i r 

2 


The program iterates at the same time step until the com- 
positions in the column have converged. The details of 
the formulation of the Langmuir isotherm appears in 
Appendix B. 

On preliminary studies done with this program, 
it was noticed that the indices , 11 and 12 *, did not 

S ee~pr o'gr* anTIT stl ng7 Appendix C 
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change during the iterations. It was decided to perform 
calculations only from II to 12 for the subsequent itera- 
tions, hence eliminating unnecessary computation time 
without losing accuracy. This simple modification resulted 
in approximately a further 10% savings of computation time. 
Thus, with both modifications, an overall reduction or 
about 42% of the computation time is achieved. 

For the n-heptane system on the Chromosorb-102 
at 175°C (Meisch, 1973), three different values of the 
Langmuir constant, K (0, 5, and 10) were studied. These 
simulations are depicted in Figures 10,. 11, a r»u 12 res- 
pectively. For this system, n-heptane on the Chromosorb- 
102 at 175°C, the input pulse was not as .arrow as one 
would like. Because of this wide input, many data points 
(over two hundred) are needed to describe it. With the 
discreteness of the input data points, an interpolation 
routine is needed to obtain input compositions at inter- 
mediate values. The cubic spline interpolation routine 
is used. With so many data points, the determination of 
the cubic spline coefficients becomes a problem. The 
subroutine, ICS1CE ( IMSL, 1973) can only handle a limited 
number of data points. When this number is exceeded, 
underflows are the result. The subroutine, ICSlCE, will 
determine the cubic spline coefficients but when they are 
used oscillations appear in the interpolations. These 
oscillations then appear in the output of the simulation 
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as can be seen by the oscillations at the peaks of the 
simulation. For smaller sets of data points, no or ‘..illations 
appear at the peaks of the simulation. However, when 
using a large set of data points to represent the input, 
another interpolation scheme should be considered. 

The value of K equal to zero corresponds to 
the linear model and the simulated results from both 
programs match identically. The three simulations of 
Figures 10, 11, and 12 can be compared with the actual 
data in Figure 6. For the Langmuir constant equal to 
zero the peak is not as sharp and there is not much 
spreading in the tail. As K is increased the rise of 
the peak is sharper and approximates the actual peak 
well. There is also more spreading in the tail. Unfor- 
tunately, the peak times of the simulations deviate 
more from the observed data peak time. Table 3 shows this 
trend. This indicates that the initial estimate of the 
thermodynamic parameter, mRo, is wrong. This suggested 
that another way of determining the parameter, mRo, had 
to have been found. This concept is discussed in the next 
section. 
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TABLE 3 


COMPARISON OF THE SIMULATION PEAK TIMES 
OF THE OUTPUT FOR VARIOUS VALUES OF THE 
LANGMUIR CONSTANT, K 
(N-HEPTANE ON CHROMOSORB-102 AT 175°C) 


K 


Peak Time 
(Dimensionless 
Time) 


0 

5 

10 


28.920 

25.400 

23.800 


ACTUAL 


29.005 







Previous research done by Lavoie (1974) and 
Meisch (1969) have shown that the parameter, mRo, is 
dependent on composition. Lavoie ran experiments with 
n-heptane at 200°C on the Chromosorb-102 by varying the 
sample size. He obtained an mRo for each experiment. 

In the current studies, it was assumed that the mRo 
corresponding to the dilute sample applied for all sample 
sizes. The Langmuir constant could then be varied until 
the peak times of the actual data and of the simulation 
match. 

A numerical scheme is needed todetermine the 
Langmuir constant. The secant method (Conte and de Boor, 
1972) is used. Coding problems had to be dealt with when 
this secant method was implemented into the program. 

These problems arose from the oscillations produced by 
the interpolation of the input data at the entrance of 
the column. Fortunately, these oscillations damped out 
well before the peak of the output occurred. To determine 
the peak times of the simulation, the pulse maximum was 
required. With the early oscillations of the simulation, 
one would not get the true maximum of the pulse so that 
a level, TEST2, was introduced into the program. Nothing 
was done below this level. As soon as the compositions 





n 


o 

0 



11 

f 

*. I 


0 

0 

0 

II 




became greater than this level, a search was made for the 
maximum peak. In this way the true peak time of the sim- 
ulation could be obtained. Figure 13 shows this general 
procedure . 

Another serious problem occurred due to these 
oscillations produced by the cubic spline interpolation 
of the input data. The pulse should always be propagating 
down the column. Consequently, it should never back track. 
In terms of the nomenclature of the program the indices, 

11 and 12, should never decrease. This was not the case 
for small values of mRo, so that a nev; feature was re- 
quired to avoid the divide checks (division by zero) that 
occurred due to this back tracking. The values of II and 

12 after a successful set of iterations are stored as 
HOLD and I20LD, respectively. During the course of the 
program, a check is constantly being made zo that the 
new valuesof II and 12 be atleast equal to their old 
values. This would eliminate any reversing of the indices 
and avoid any computer divide checks. 

In Figure 14 an organization chart is presented 
to show the input and output of the Langmuir program. 

A second n-heptane system was studied. This is 
one of the sample size mRo dependence experiments run by 
Lavoie. The actual chromatograph for n-heptane on the 
Chromosorb-102 column operating at 200°C is shown in 
Figure 15 and the equilibrium adsorption model chromatograph 
is presented in Figure 16. The equilibrium adsorption 
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Time 


Figure 13 Determination of Peak Time for 
Langmuir Program 








Figure 14 Organization of Langmuir Isotherm Program 







Figure 16 Chromatograph for n-IIeptane on Chromosorb-102 
Using the Equilibrium Adsorption Model 




model curve is very sharp and the peak is about three times 
as great as the actual data. The simulation curve exhibits 
almost no spreading as opposed to the high degree of tail- 
ing in the experimental data. The compositions for this 
system are much larger than encountered before and the 
equilibrium adsorption model poorly predicts the output 
as expected. 

For the simulation using the Langmuir isotherm, 
the parameter, mRo, correspoding to the sample size of, 

0.2 microliters was used to simulate data for the system 
having the larger sample size of 3.0 microliters. The 
Langmuir simulation for n-heptane on the Chromosorb-102 
at 200®C with an mRo » 0.0482 and the Langmuir constant 
of 4.8 is shown in Figure 17. This simulation approxi- 
mates the actual data quite well. It rises quite sharply 
like the actual output peak and exhibits similar tailing. 
The peak height is still slightly higher than the actual 
data but nowhere near as great as the linear isotherm 
prediction. To further improve the peak height, a slightly 
smaller mRo would seem appropriate. A value of mRo equal 
to 0.0382 (K»12.2) produced a simulation with a slightly 
smaller maximum. 'This simulation chromatograph is 
shown in Figure 18. A third and final simulation was made 
with an intermediate value of mRo equal to 0.0426 (K»8.6) 
and its simulation is presented in Figure 19. These three 
simulation peak heights are compared in Table 4. 

















TABLE 4 


COMPARISON OF THE SIMULATION PEAK HEIGHTS FOR 
VARIOUS MRO'S WITH THE ACTUAL I EAK HEIGHT 
(N-HEPTANE ON CHROMOSORB-102 AT 200°C) 


Thermodynamic 

Langmuir 

Peak Height 

Parameter, mRo 

Constant, K 

(mole fraction) 

0.0382 

12.2 

0.0308 

0.0426 

8.6 

0.0333 

0.0482 

• 

OO 

0.0430 


ACTUAL DATA 


0.0353 


With successful simulation of the n-heptane 
system on the Chromo sorb-10 2 column at 200°C, attention 
was directed to the n-heptane system on the Chromosorb-102 
column at 17 5° C. Since there are no data for a dilute 
sample, an estimate of the parameter, mRo, was made. With 
this approximate value of the parameter, difficulties 
in convergence of the Langmuir constant v/ere encountered. 
The oscillations from interpolation are causing the 
convergence problems of K. No convergence was obtained 
for this system. These problems were a result of the 
problems of the cubic spline interpolation. Another 
interpolation scheme should be investigated. A different 
numerical scheme might be studied which would converge 
faster to the appropriate Langmuir constant. 
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This investigation has been conducted in con- 
junction with the group effort to define fundamental design 
criteria necessary for an optimal design of a combination 
gas chromatograph-mass spectrometer. Specifically, this 
report has dealt with the formulation of a modification 
of the equilibrium adsorption model with the incorporation 
of the Langmuir isotherm and subsequent efforts to ascertain 
the merits of the numerical technique known as Crank-Nicol- 
son as a technique worthy of use in the simulation of a 
nonlinear gas chromatograph model. 

Previous work suggested the formulation of a 
model which took into account the variation of the thermo- 
dynamic parameter, mRo, with composition. The Langmuir 
isotherm has been used to account for the fact that there 
is a saturation of adsorbent sites at higher concentrations. 
The analysis of the equilibrium adsorption model using 
the Langmuir isotherm indicates the characteristics of the 
actual data are more adequately predicted than with previous 
models . 

The finite difference method of Crank -Nicolson 
has shown to have less instabilities than the orthogonal 
collocation method studied earlier. The Crank-Nicolson 
method can be modified to shorten computation time with 
no significant loss in precision, with the two modifications 
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implemented into the algorithm the computation time was 
reduced approximately 42%. 

The interpolation method used on the discrete 
set of input points exhibited the property that for a 
large set of data points (a wide pulse) underflows in 
calculating the cubic spline coefficients resulted. These 
underflows manifested themselves later in the simulation 
by producing oscillations at the peak of the output. For 
a case with few data points no oscillations occurred. 
Therefore, an attempt should be made in the future to find 
a better interpolation scheme. 

The incorporation of the secant method into 
the simulation model to match the peak times of the simu- 
lation for a given mRo with the actual peak time resulted 
in good approximations to the actual data. For a better 
approximation to the tailing and the steep rise of the 
peak an mRo corresponding to a dilute sample was used. 

The parameter, mRo, could be further varied to better fit 
the actual data. Future work could incorporate another 
secant method or least squares fit to vary mRo which in 
turn would determine the Langmuir constant, K, to match 
the peak heights and peak occurrences of the simulation 
with those of the actual data. This results in what amounts 
to as a two constant curve fitting scheme. 

When the numerical scheme was applied to the 


first system studied there was no convergence for the 
Langmuir constant. This problem was the result of the 
oscillations produced from the interpolation of the input 
using the cubic spline coefficients. This suggests again 
the importance of finding a better numerical scheme for 
interpolating the input data. A better convergent numer- 
ical technique not as dependent on the oscillatons of the 
peak should be investigated. 

The Crank-Nicolson method solves the nonlinear 
equation quite well but even with the modifications made 
to it an excessive amount of computation time is needed. 
Depending on the system, between 25 and 60 minutes of 
computation time is needed. Another area of future work 
could be to go back to the method of orthogonal collocation 
and try . new trial functions and new orthogonal functions 
which would result in shorter simulations. Other numerical 
techniques might be studied also. 
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